Activity: Spatially Continuous Data III

Activity: Spatially Continuous Data III

Practice questions

Answer the following questions:

1. What is a correlogram?
2. What is the relationship between the autocovariance and the semivariance?
3. Describe the elements of a semivariogram.
4. Why is it important to consider the number of pairs used in the calculation of the semivariance?

Learning objectives

In this activity, you will:

1. Calculate and plot empirical semivariograms.
2. Estimate and plot theoretical semivariograms.
3. Discuss the results of variographic analysis.

Suggested reading

• Bailey TC and Gatrell AC [-@Bailey1995] Interactive Spatial Data Analysis, Chapters 5 and 6. Longman: Essex.
• Bivand RS, Pebesma E, and Gomez-Rubio V [-@Bivand2008] Applied Spatial Data Analysis with R, Chapter 8. Springer: New York.
• Brunsdon C and Comber L [-@Brunsdon2015R] An Introduction to R for Spatial Analysis and Mapping, Chapter 6, Sections 6.7 and 6.8. Sage: Los Angeles.
• Isaaks EH and Srivastava RM [-@Isaaks1989applied] An Introduction to Applied Geostatistics, CHAPTERS. Oxford University Press: Oxford.
• O’Sullivan D and Unwin D [-@Osullivan2010] Geographic Information Analysis, 2nd Edition, Chapters 9 and 10. John Wiley & Sons: New Jersey.

Preliminaries

For this activity you will need the following:

• An R markdown notebook version of this document (the source file).

• A package called geog4ga3.

It is good practice to clear the working space to make sure that you do not have extraneous items there when you begin your work. The command in R to clear the workspace is rm (for “remove”), followed by a list of items to be removed. To clear the workspace from all objects, do the following:

rr rm(list = ls())

Note that ls() lists all objects currently on the worspace.

Load the libraries you will use in this activity (load other packages as appropriate).

library(tidyverse)
library(gstat)

package 㤼㸱gstat㤼㸲 was built under R version 3.5.3
Attaching package: 㤼㸱gstat㤼㸲

The following object is masked from 㤼㸱package:spatstat㤼㸲:

    idw

library(spdep)
library(geog4ga3)

Load dataset:

data("aquifer")

Convert to a SpatialPointsDataFrame:

aquifer.sp <- aquifer
coordinates(aquifer.sp) <- ~X+Y

The data is a set of piezometric head (watertable pressure) observations of the Wolfcamp Aquifer in Texas (https://en.wikipedia.org/wiki/Hydraulic_head). Measures of pressure can be used to infer the flow of underground water, since water flows from high to low pressure areas.

These data were collected to evaluate potential flow of contamination related to a high level toxic waste repository in Texas. The Deaf Smith county site in Texas was one of three potential sites proposed for this repository. Beneath the site is a deep brine aquifer known as the Wolfcamp aquifer that may serve as a conduit of contamination leaking from the repository.

The data set consists of 85 georeferenced measurements of piezometric head. Possible applications of interpolation are to determine sites at risk and to quantify uncertainty of the interpolated surface, to evaluate the best locations for monitoring stations.

Activity

1. Obtain and plot the empirical semivariogram for the head in the Walker Lake dataset. How would you interpret this semivariogram?

variogram_z <- variogram(H ~ 1, data = aquifer.sp)
ggplot(data = variogram_z, aes(x = dist, y = gamma)) +
  geom_point() + 
  geom_text(aes(label = np), nudge_y = -1500) +
  xlab("Distance") + ylab("Semivariance")

2. Estimate trend surface of your choice, and obtain and plot an empirical semivariogram using the residuals. How would you interpret this semivariogram?

trend.l <- lm(formula = H ~ X + Y, data = aquifer)
summary(trend.l)

Call:
lm(formula = H ~ X + Y, data = aquifer)

Residuals:
    Min      1Q  Median      3Q     Max 
-367.00 -161.41  -30.74  148.23  651.17 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2591.3266    38.9636   66.51   <2e-16 ***
X             -6.7519     0.3439  -19.64   <2e-16 ***
Y             -5.9862     0.4066  -14.72   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 203.3 on 82 degrees of freedom
Multiple R-squared:  0.892, Adjusted R-squared:  0.8894 
F-statistic: 338.8 on 2 and 82 DF,  p-value: < 2.2e-16

u.p <- seq(from = 0.0, to = 1.0, by = 0.05)
v.p <- seq(from = 0.0, to = 1.0, by = 0.05)
preds <- predict(trend.l, newdata = aquifer, se.fit = TRUE, interval = "prediction", level = 0.95)
z.p <- matrix(data = preds$fit[,1], nrow = 21, ncol = 21, byrow = TRUE)

data length [85] is not a sub-multiple or multiple of the number of rows [21]

z.p_l <- matrix(data = preds$fit[,2], nrow = 21, ncol = 21, byrow = TRUE)

data length [85] is not a sub-multiple or multiple of the number of rows [21]

z.p_u <- matrix(data = preds$fit[,3], nrow = 21, ncol = 21, byrow = TRUE)

data length [85] is not a sub-multiple or multiple of the number of rows [21]

trend.plot <- plot_ly(x = ~u.p, y = ~v.p, z = ~z.p, 
        type = "surface", colors = "YlOrRd") %>% 
  add_surface(x = ~u.p, y = ~v.p, z = ~z.p_l, 
              opacity = 0.5, showscale = FALSE) %>%
  add_surface(x = ~u.p, y = ~v.p, z = ~z.p_u, 
              opacity = 0.5, showscale = FALSE) %>% 
  layout(scene = list(
    aspectmode = "manual", aspectratio = list(x = 1, y = 1, z = 1)))
trend.plot #%>% add_markers(data = df, x = ~x, y = ~y, z = ~z, color = ~residual)

3. Estimate and plot a theoretical semivariogram model for the residual variogram.

aquifer$residual1 <- ifelse(trend.l$residuals > 0, "Positive", "Negative")
ggplot(data = aquifer, 
       aes(x = X, y = Y, color = residual1)) +
  geom_point() +
  coord_equal()

4. How would you use the information provided by the variographic analysis above to improve your predictions of the field?